Event-by-event study of prompt neutrons from 239 Pu(n, /) 
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Employing a recently developed Monte-Carlo model, we study the fission of 240 Pu induced by 
neutrons with energies from thermal to just below the threshold for second chance fission. Current 
measurements of the mean number of prompt neutrons emitted in fission, together with less accurate 
measurements of the neutron energy spectra, place remarkably fine constraints on predictions of 
microscopic calculations. In particular, the total excitation energy of the nascent fragments must be 
specified to within 1 MeV to avoid disagreement with measurements of the mean neutron multiplicity. 
The combination of the Monte-Carlo fission model with a statistical likelihood analysis also presents 
a powerful tool for the evaluation of fission neutron data. Of particular importance is the the fission 
spectrum, which plays a key role in determining reactor criticality. We show that our approach can 
be used to develop an estimate of the fission spectrum with uncertainties several times smaller than 
current experimental uncertainties for outgoing neutron energies of less than 2 MeV. 
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I. INTRODUCTION 

The quest for a fundamental theory of fission began 
with the 1939 seminal work of Bohr and Wheeler [l|, 
the same year as this phenomenon was discovered by 
Hahn and Strassmann [2[ and interpreted by Meitner and 
Frisch [3] . Bohr and Wheeler used the liquid-drop model 
to make predictions that were remarkably realistic given 
the paucity of available data. The current theoretical de- 
scriptions of fission reflect the complexity and richness 
revealed over 70 years of experimental studies, empha- 
sizing the multi-dimensional, dynamic, and microscopic 
aspects. In particular, a refined version of the liquid drop 
model that includes a finite interaction range and quan- 
tum shell corrections has formed the basis for extensive 
calculations of the potential-energy surfaces associated 
with the multidimensional shape of fissioning nuclei (see 
Refs. [1, @ and references therein). Concurrently, a pro- 
gram is underway to develop a fully microscopic treat- 
ment of fission in terms of a quantum many-body treat- 
ment of protons and neutrons subject to an adjustable 
effective (in- medium) interaction [1, 0, Hi- 

Despite the many theoretical advances, there is not 
yet a quantitative theory of fission. This is unfortunate 
because nuclear fission remains important to society at 
large due to its many practical applications, including en- 
ergy production and security. For example, reactors and 
other critical systems demand that neutron growth be 
known to about the 0.1% level for model simulations to 
be reliable. In such cases, scattering experiments are in- 
sufficiently accurate, requiring reliance on more inclusive, 
higher statistics integral critical assembly experiments. 

Furthermore, in the last few years efforts have been un- 
derway to develop systems capable of detecting concealed 
nuclear material. These applications place entirely differ- 
ent demands on fission models by attempting to exploit 
specific information carried by particles resulting from 
fission. Thus there is a need for a fission description that 
accounts for particle correlations and fluctuations on an 



event-by-event level. Such a description, employing a 
model incorporating the relevant physics with a few key 
parameters, compared to the pertinent data through a 
statistical analysis, presents a potentially powerful tool 
for bridging the gap between current microscopic mod- 
els and important fission observables and for improving 
estimates of the relatively gross fission characteristics im- 
portant for applications. This type of approach also pro- 
vides a means of using readily measured observables to 
constrain our understanding of the microscopic details of 
fission. 

Relatively recently, Lemaire et al. [9j implemented a 
Monte-Carlo simulation of fission fragment statistical de- 
cay by sequential neutron emission for spontaneous fis- 
sion of 252 Cf and thermal fission of 235 U. That work 
demonstrated how fission event simulations, in conjunc- 
tion with experimental data on fission neutrons and 
physics models of fission and neutron emission, can be 
used to predict the neutron spectrum and to validate 
and improve the underlying physics models. 

In the present work, we have implemented a concep- 
tually similar approach and applied it to calculate the 
sequential neutron emission for the neutron induced fis- 
sion of 240 Pu. Specifically, we have adapted the recently 
developed fission event generation model FREYA [29( to 
calculate the production and decay of fission fragments 
and used maximum-likelihood analysis to estimate prop- 
erties of the emitted fission neutrons and their correlation 
coefficients. To our knowledge, such correlations have not 
been extracted before for fission neutrons in a physics- 
based Monte-Carlo simulation. The detailed statistical 
analysis presented here is essential for developing a more 
quantitative understanding of fission and obtaining bet- 
ter evaluations of fission data for various applications. 

First, in Sect. HQ we present the framework for the 
statistical analysis employed for obtaining estimates of 
the model parameters and the neutron observables, as 
well as the correlations between the various quantities of 
interest. We then discuss in Sect. IIIII the experimental 
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data used in this work with a particular emphasis on 
experimental uncertainties. Subsequently, in Sect. lIVl we 
describe the physics ingredients of the FREYA simulations. 
Finally, in Sect. [V] we present calculated results for the 
239 Pu(?i, /) neutron spectrum and other observables for 
incident neutron energies, E n , from 0.5 to 5.5 MeV. 



II. STATISTICAL METHOD 

Here we briefly describe the statistical method used for 
determining model parameters and reaction observables. 

There are a number of different techniques for estimat- 
ing model parameter values and although their relative 
merits are being vigorously debated they often differ very 
little in their actual results. Our present analysis is in- 
spired by the general inverse problem theory developed 
by Tarantola [301 ] . 

We introduce a number of model parameters {a k } (de- 
fined in Sect. IIV[) . Since the theory does not, a priori, 
specify the parameter values, we assume that the param- 
eter values are uniformly distributed over a reasonable 
interval in parameter space. For a specified set of pa- 
rameter values {ai }j we generate a large sample of 
fission events from which we then extract the particu- 
lar observables of interest, {Ci}. These calculated values 
are then compared with the corresponding experimental 
values, {£i}. 

Specifically, for each parameter set {ot^ } we calculate 
the x 2 deviation of the calculated observables from their 
measured values, 
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Here {&{} are the uncertainties in the experimental values 
Division by these quantities ensures that well-measured 
observables carry more weight than those that are poorly 
measured. 

The key feature of the method [3(| is that a likelihood 
is assigned to each particular set m of model parame- 
ter values based on how well the corresponding model 
calculation reproduces the experimental results, 



r ( m ) i 
w{a k } oc 



(2) 



This quantity is then taken as the relative probability 
that those parameter values are the "correct" ones. In 
this manner, one may define a probability density in the 
space of model parameters, P{a k } = w{a k }/W, where 

W = J2 ln w ™ 1S t ne sum °f an the weights. 

Once the probability density of model parameter val- 
ues has been obtained, their corresponding statistical dis- 
tribution of the observables can readily be calculated. 
Thus the best estimate for the model parameter values, 
{(5;^}, is given by the likelihood-weighted average, 
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The last relation indicates that the best estimate is ap- 
proximately equal to the most likely value i.e. the 
value having the largest likelihood. The covariances 
among the parameter values are similiarily calculated, 
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The diagonal elements, a kk — &\, are the usual vari- 
ances with a k the standard deviations of the parame- 
ter values and represent the squares of the uncertain- 
ties on the values of the individual model parameter 
a k - The off-diagonal elements give the covariances be- 
tween two model parameters. It is often more instruc- 
tive to employ the associated correlation coefficients, 
Cfcfc' = o kk > /[d k d k >]. 

An analogous procedure can be carried out to obtain 
best estimates for the various calculated quantities, i.e. 
for the observables {d}. Thus, if C\ m) = C t {a[ m) } de- 
notes the value of Ci calculated with the particular pa- 
rameter values {a k ™ }, then the best estimate for the 
observable Ci is given by 
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The last relation expresses the fact that the best estimate 
is approximately equal to the most likely result, i.e. the 
result obtained with the most likely parameter values. 

Covariances between different observables, {Ci}, are 
calculated as 



= -co (c, -<*•)>- 
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The diagonal elements are the squares of the standard 
deviations, {&i\, of the calculated values {Ci} resulting 
from uncertainties in the model parameter values. Here 
Cij = aij/[aidj\ are the correlation coefficients between 
the observables Ci and C 7 - . 
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FIG. 1: (Color online) The evaluated ENDF/B-VII data for 
the average prompt neutron multiplicity 77 as a function of the 
incoming neutron e nergy E n , together with the experimental 
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Data Set 


AT 


£ min (MeV) 


£ max (MeV) 


a (MeV) 


b (MeV) 


X 2 /N 


Abramson [21] 


95 


0.55 


14.253 


1.042 


0.5294 


2.073 


Aleksandrova [20] 


54 


1.503 


11.128 


0.914 


0.5033 


13.474 


Aleksandrova [20] 


19 


1.5 


11 


0.917 


0.5033 


14.666 


Belov [18] 


18 


0.3 


7 


0.991 


0.5033 


0.868 


Conde [17] 


13 


0.3 


7.5 


0.975 


0.5365 


1.121 


Knitter [19] 


183 


0.28 


13.87 


1.030 


0.5040 


1.529 


Nefedov [22] 


65 


0.139 


7.15 


1.023 


0.5053 


0.765 


Starostov [25] 


65 


3.007 


11.2 


0.995 


0.5288 


3.890 


Werle [24] 


79 


0.104 


9.5 


1.035 


0.5263 


4.244 


Staples (0.5 MeV) [23] 


68 


0.615 


16 


1.026 


0.5005 


4.067 


Staples (1.5 MeV) [23] 


59 


1.7 


15.2 


1.009 


0.5025 


8.137 


Staples (2.5 MeV) [23] 


51 


2.77 


14.4 


1.0276 


0.5025 


4.018 


Staples (3.5 MeV) [23] 


38 


4.07 


13.8 


1.0354 


0.5025 


8.033 



TABLE I: For each data set is listed the number of points N, the minimum and maximum measured outgoing neutron energies, 
the fitted Watt parameters a and 6, and the associated x 2 P er degree of freedom. 



In principle, the best estimate for the observables {Ci} 
is neither that resulting from using the most likely pa- 
rameter values {o^} nor that calculated with the best 
estimate of the model parameters, In our appli- 

cations the distinction between the different estimates is 
mostly one of principle since the different estimates yield 
practically identical results. We shall generally adopt the 
observable values calculated with the optimal parameter 
values, {a 9 }, as our estimate while the associated uncer- 
tainties and correlations will be obtained on the basis of 
the entire ensemble, as expressed in Eq. ([5]). 



III. EXPERIMENTAL DATA 

We discuss here the experimental data used in our 
study. 



A. Mean neutron multiplicity 



The mean number of prompt neutrons emitted follow- 
ing neutron- induced fission of 239 Pu has been measured 
in a number of experimen ts JTol . [Til H3 , [T3L ITij and was 
reviewed by Fort et al. jT^j • Figure Q] shows a selec- 
tion of this data as well as the associated ENDF/B-VII 
evaluation [l6|]. We employ the ENDF evaluation as an 
approximate average of the experimental numbers and 
we assign a 0.5% uncertainty to V. 



B. Prompt neutron spectrum 

Our statistical analysis will also incorporate the mea- 
sured prompt neutron spectrum p^[l^ .[lJ.[20l.[2lT[2^.[23] 
as given in the EXFOR/CSISRS database. Wherever the 
experimental uncertainties are not given we have used an 
uncertainty of 5% in the calculation of x 2 - This is likely 
an under- estimate of the real uncertainty. The various 
data sets are shown in Fig. [21 The bulk of the data are ob- 
tained for low incident neutron energies, E n < 0.5 MeV. 
The remaining data have been taken by Staples et al. [23j 
at E n = 0.5, 1.5, 2.5, and 3.5 MeV. 

The data in the top panel of Fig. [5] were taken for 
incident energies below 0.5 MeV and are not absolutely 
normalized. In order to compare the data sets with each 



other and with our calculated spectra, we normalize all 
data sets to unity (while preserving the spectral shapes). 
For this purpose, we fit the observed energy spectra to a 
Watt spectrum, 

^ = N e- E ' a smhy/2E]b, (7) 

where the normalization Nq is determined by demanding 
that the integral over E yield unity. Table U lists for each 
data set the number of data points, the minimum and 
maximum neutron energies observed, the Watt parame- 
ters a and b obtained by the fitting procedure, and the 
associated minimum x 2 per degree of freedom. The value 
of a is 1 MeV within the uncertainties of the fits for 
all but the Aleksandrova sets where a sa 0.91 MeV. The 
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value of b is 0.50-0.56 MeV in all cases. 

The data on the neutron spectra cover a wide energy 
range, 0.1 < E < 14 MeV. In the lowest E range, E < 0.5 
MeV, the neutron yields generally increase with E, reach- 
ing a maximum somewhere between 0.5 and 1 MeV, and 
decreasing again above 1 MeV. There is significant dis- 
agreement between the data sets in this energy region. In 
particular, the data of Belov et al. [l8j . Werle et al. [24| . 
and Abramson et al. (2lj have relatively large uncertain- 
ties and include points noticeably higher than the re- 
maining data. Curiously, the peak of the E n — 0.5 MeV 
spectrum from Staples et al. [23| is significantly narrower 
than those of the other data sets. At higher outgoing 
energies, E > 2 MeV, all the data sets closely follow 
each other, except for those from Aleksandrova et al. [2(| 
which are systematically lower. Indeed, the Aleksandrova 
sets are rather poorly represented by the Watt fits, hav- 
ing the largest x 2 per degree of freedom, see Table |U 

Some of the discrepancies between the data sets may be 
due to the incompleteness of the individual sets in parts 
of the energy range. For example, the Aleksandrova [2fj| 
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FIG. 2: (Color online) The measured prompt neutron energy 
spectra, normalized to unity, as a function of outgoing neutron 
energy for low incident energies from Refs. [13, GIL Gil, H3, Hi], 
I22I I23I. I24I |25|] (upp er p anel) and for a wider range of incident 
energies from Ref. [23j ] (lower panel). 
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FIG. 3: (Color online) The average TKE as a function of the 
heavy fragment mass Ah, from Refs. [26l. [27t [28(] . 



and Starostov [25j sets are only available for E above 1.5 
and 3 MeV, respectively, so that the Watt fits may match 
the high energy tail of the spectrum but cannot represent 
the peak region and below. Similarly, sets that cover 
the region E < 7 MeV may not give as good fits to the 
high-energy part of the spectrum. When E n > 1.5 MeV, 
the minimum outgoing energy E measured by Staples et 
al. [23j (shown in the lower panel of Fig. ^ is always 
greater than E n . Thus these data sets do not provide 
much information on the softer part of the spectrum and 
the back extrapolation by means of the Watt form is 
somewhat unreliable since the measured hard spectra do 
in fact not fit a Watt shape very well, as reflected by the 
large values of \ 2 in Table Q] 



C. Fission fragment energies 

Several measurements of the total kinetic energy 
(TKE) of the two fission fragments can be found in the 
literature. Figure [3] shows the principal measurements of 
the mean TKE as a function of the mass number of the 
heavy fragment, Ah, which were made by Wagemans et 
al. [ll], Nishio et al. (27[, and Tsuchiya et al. [28]. (The 
mass number of the heavy fragment is found by simul- 
taneously measuring the velocities and energies of both 
fragments [27| . No experimental uncertainties are given 
for these results, neither for the mass number nor for the 
reported TKE.) The data exhibit a significant dip near 
symmetry and fall off steadily for large asymmetries, re- 
sulting in a maximum at Ah ~ 133. The different data 
sets generally agree well for large Ah but they exhibit a 
significant spread near symmetry. Furthermore, Ref. [27tj 
also provides the full-width at half-maximum (FWHM) 
of the TKE distribution at selected values of Ah ■ These 
also decrease at large Ah, reflecting the fact that the 
TKE spectrum softens, presumably because the mutual 
Coulomb repulsion between the two nascent fragments 
decreases with larger asymmetry. 
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IV. GENERATION OF FISSION EVENTS 

We have adapted the recently developed fission model 
FREYA [2!| for the present purpose of calculating the neu- 
tron spectrum in terms of a set of well-defined model pa- 
rameters. Since this is the first practical application of 
FREYA, we describe its main physics ingredients below. 

The code follows the temporal sequence of individual 
fission events from the initial excited fissionable nucleus, 
240p u * m p re g en ^ cage, through a scission configura- 
tion of the two nascent fragments, to the subsequent neu- 
tron evaporation from the fully accelerated fragments. 
The competition between fission and neutron emission 
from the fissioning nucleus (2 nd chance fission) has not 
yet been implemented in the code. Consequently, we re- 
strict our discussion to energies below 5.5 MeV. 



A. Fission mass and charge partition 

The fission process is initiated when a neutron with a 
specified initial energy E n is absorbed by a fissile nucleus 
to form a compound nucleus A Z with a certain excitation 
energy. The compound nucleus subsequently splits into 
a heavy fragment Ah Zh and a complementary light frag- 
ment Al Zl. In its present early form, FREYA selects the 
mass and charge partitions on the basis of existing ex- 
perimental data. For the present study, we use fits to the 
thermal and fast 239 Pu(n, /) fission product mass yields 
measured by England and Rider [3lT ] in combination with 
the charge distributions obtained by Reisdorf et at [HJ . 

The fits assume that the mass product yields Y(A p ) ex- 
hibit three distinct fission modes that can be represented 
in terms of suitable gaussians, 



Y(A P ) = S 1 (A P ) + S 2 (A P ) + S L (A P ) 



(8) 



The first two terms result from asymmetric fission modes 
associated with the spherical shell closure at N = 82 and 
the deformed shell closure at N = 88 respectively, while 
the last term results from a symmetric, so-called super- 
long, mode which is relatively insignificant 33] . The spe- 
cific forms of these terms are 



The values of the parameters in the fits to Y{A p ) are 
given in Table [TT| for either thermal or fast fission. The 
normalization is chosen such that ^2 A Y(A) = 2 since 
each event leads to two products. Consequently we have 
2Ni + 2N2 + Nl = 2, apart from a negligible correction 
because A p is discrete quantity bounded both from be- 
low and above. It should be noted that the symmetric 
component contributes only 1-2 per mille of the yield. 

While these fits are to the fission product yields, the 
FREYA simulation requires fission fragment yields, i.e. the 
probability distribution for obtaining a given mass par- 
tition at scission, before neutron evaporation has begun, 
We take A w \Aq , but keep the displacements D, t and 
the widths cr* unchanged. We use the thermal fits for 
E n < 1 MeV and the fast fits for 1 < E n < 5.5 MeV, the 
highest value of E n considered here. The change in the 
fit parameter values with incident neutron energy should, 
of course, be more continuous than we have implemented 
here but the change is most important in areas where the 
yields are low: the tails of the gaussians where fission is 
most asymmetric and in the case of symmetric fission. 
At even higher energies, symmetric fission (the Sl com- 
ponent) grows increasingly important, filling in the dip 
at symmetry. The width cr 2 also increases, broadening 
the asymmetric tails. 

The resulting fits are compared to the data in Figs. Q] 
and [5J The agreement with the tabulated percentage 
yields is quite good, especially in the regions where the 
yields are highest and which thus contribute the greatest 
number of events. Equation © does not perfectly de- 
scribe the tails at high and low fragment mass. We have 
also tried a fit with 5 independent gaussians, e.g. allow- 
ing Ni, Di and <jj to vary independently on the low and 
high sides of A, and found that the fit does not signifi- 
cantly improve as a result. We note also that the width 
of the Sl component is not as large as found in other 
actinides where the yields have been decomposed in a 
similar fashion (34|. 

Once the gaussian fit has been fixed, it is straightfor- 
ward to make a statistical selection of the fragment mass 
number Af. The mass number of the partner fragment is 
then readily determined since we assume Al + Ah = Aq. 
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V 277(7,; 
for i = 1,2 and 



-(A-A-Di) 2 /2of 1 e -(A-A+Di) 2 /2<r? 



(9) 



2-kgl 



(10) 



Here A = \{A 



where Ao = 240 is the mass num- 
ber of the fissioning nucleus and v is the average total 
multiplicity of evaporated neutrons. (While there ex- 
ist more detailed data for e.g. 235 U(rt, /) that give the 
yields as a function of both mass and total kinetic energy, 
Y (A p , TKE), for several values of E n (34[, such data are 
not yet available for Pu.) 



Parameter 


Thermal 


Fast 


A 


118.5 


117.5 


Ni 


0.7574 


0.7355 


Dj 


20.81 


20.96 


<Tl 


5.626 


5.711 


N 2 


0.2417 


0.2623 


D 2 


14.95 


15.14 


02 


2.546 


2.627 


N L 


0.0018 


0.0044 




1.824 


2.511 



TABLE II: The fit parameters of the three fission modes for 
thermal and fast neutron-induced fission on 239 Pu. 
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240 Pu*, and ground-state masses of the two fragments, 

Qlh = M ( 240 Pu*) - M L - M H . (12) 

FREYA takes the required nuclear ground-state masses 
from the compilation by Audi et al. J35|, supplemented 
by the calculated masses of Moller et al. [3(| where no 
data are available. The Qlh value for the selected fis- 
sion channel is then divided up between the total kinetic 
energy (TKE) and the total excitation energy (TXE) of 
the two fragments. The specific procedure employed is 
described below. 

First, the average value of TKE is determined on the 
basis of the Coulomb potential between the two frag- 
ments at scission, 



FIG. 4: (Color online) The fission product yield as a func- 
tion of fragment mass for thermal fission. The data are from 
Ref. [3l| ] while the curves are a five-gaussian fit to the data. 
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FIG. 5: (Color online) Same as Fig. [4] but for fast-neutron 
induced fission. 



The fragment charge, Zf, is selected subsequently. For 
this we follow Ref. [32j and employ a gaussian form, 



P Al {Z f ) 



oc e 



-(Zf-Z f (A f )) 2 /2a 2 z 



(11) 



with the condition that \Zf — Zf(Af)\ < 5az- The cen- 
troid is determined by requiring that the fragments have, 
on average, the same charge-to-mass ratio as the fission- 
ing nucleus, Zf(Af) = AfZ^/A^. The dispersion is the 
measured value, oz = 0.5 [32|]. The charge of the com- 
plementary fragment then follows using Zl + Zh = Zq. 



B. Fragment energies 

Once the partition of the total mass and charge among 
the two fragments has been determined, the Q value asso- 
ciated with that particular channel follows as the differ- 
ence between the mass of the excited compound nucleus, 



TKE 



Zr.Z 



CL + CH 



d-LH 



(13) 



In the scission configuration, the two nascent fragments 
are assumed to have spheroidal shapes and be positioned 
coaxially with a tip separation of cLlr- The associated 



major axes are c, = roA^ 3 /[l — ^e(Zi, Aj)] with r 
fm. We use the values for the spheroidal deformation 
parameter e(Zi, Aj) calculated in Ref. [36| which include 
shell effects. The denominator of Eq. (fl3|) is thus the 
distance between the centers of the two fragments and 
the above expression represents the monopole-monopole 
term of the mutual Coulomb interaction energy. 

The tip separations {cIlh} are important parameters 
in the model since they determine the (average) fragment 
kinetic energies and hence, by energy conservation, also 
the total fragment excitation that is available for neutron 
emission. Thus the neutron emission is quite sensitive 
to the specified values of {<1lh} and they deserve care- 
ful consideration. Furthermore, since the TKE is closely 
related to the Coulomb potential at scission, these pa- 
rameters contain valuable information about the scission 
configurations. 

Figure [S] shows the mean total fragment kinetic energy 
as a function of mass number of the heavy fragment as ob- 
tained by using a common tip separation dp for all fission 
channels. A comparison to the data [26|, [U [28| shows 
significant discrepancies near symmetry where the cal- 
culated TKE exhibit an enhancement whereas the data 
have a dip. 

To account for the dependence of the tip separation on 
the mass partition, we took the average of the data sets 
shown in Fig. [S] and extracted the average tip separations 
d LH shown in Fig. [3 assuming that the two fragments 
have the same charge-to-mass ratio. Near symmetric fis- 
sion, d L H is large, 7-8 fm at Ah = 120, with a steep 
drop to less than 4 fm for Ah > 132. Near symme- 
try, the plutonium fission fragments are mid-shell nuclei 
subject to strong deformation. Thus the scission config- 
uration will contain significant deformation energy and a 
correspondingly large distance between centers, resulting 
in low TKE. At Ah — 132, the heavy fragment is close 
to the doubly-magic closed shell with Zh = 50, Njj = 82 
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Heavy fragment mass number A H 

FIG. 6: (Color online) The measured average TKE as a func- 
tion of the mass number of the heavy fragment [26l . [27l . [28| 
compared to FREYA calculations with a constant tip separa- 
tion of do=4.05 fm and the average distance extracted from 

Fig.m 



and is resistant to distortions away from its spherical 
shape. However, the complementary light fragment is 
far from a closed shell and is significantly deformed, hav- 
ing thus a large value of cl which then results in a small 
tip separation d LH and a large TKE. The passage of the 
heavy fragment mass through the doubly-magic region 
results in the dip in calculated TKE around Ah ~ 130, 
see Fig. [6] 

The TKE values shown in Fig. [S] were obtained in 
experiments with incident neutrons of very low energy 
and there are no other higher-energy data to show how 
TKE(A# ) evolves with incident neutron energy. At each 
higher incident energy E n > .©thermal, we use tip separa- 
tions obtained by scaling those fitted at thermal energies, 

dLH(En) = s(E n )d LH (E t h Cima \) , (14) 

and use the common scaling factor s(E n ) as one of the 
adjustable model parameters in our fits to the neutron 
spectra. The average neutron multiplicity is very sensi- 
tive to this scale factor which, as we shall show, is greater 
than but very close to unity for the entire energy range 
studied. 

As shown in Fig. [51 the scaled tip separations lead to 
a very good agreement with the TKE data. With this 
means of fixing cIlm, the TKE is no longer overestimated 
near symmetry, leading to a better approximation of the 
individual fragment kinetic energy as well as the neu- 
tron multiplicity as a function of fragment mass, overes- 
timated and underestimated respectively with a constant 
value of d,LH: as shown in Figs. [5] and [51 The variable 
d-LH also correctly produces the dip in the single frag- 
ment kinetic energy shown in Fig. [5] The small dips in 
the fragment kinetic energy at A = 110 and 130 corre- 
spond to the dip at Ah ~ 130 in Fig. [Bj 

The overestimate of the total fragment kinetic energy 
with a constant cIlh leaves insufficient excitation energy 
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FIG. 7: (Color online) The tip separation <1lh fitted to 
the TKE values measured at thermal energies [26l . I27L l2gj ]. 
with the deformation radii extracted from the mass model of 
Ref. HI]. 

available for neutron evaporation near symmetry, result- 
ing in the near absence of neutron emission in Fig. [5] in 
this case. On the other hand, with <1lh from Fig. 
there is a peak in the neutron emission near symmetry, 
followed by a drop for A > 120, resulting in the charac- 
teristic sawtooth shape of v{A). The decrease in KE for 
these values of A gives small peaks in the neutron mul- 
tiplicity at the same values of A. Interestingly enough, 
the calculations with both fixed and variable dhH , give 
the same V even though T>{Af) is very different in the 
two cases. It is easy to see why this is true by looking at 
Figs. [4] and [9] together. Symmetric fission does not con- 
tribute significantly to the total yield, Y(Af). Most of 
the fragment yield is around Al ~ 100, Ah ~ 140. The 
variable d^n gives more neutrons for symmetric fission 
and in regions of high Ah (low Al) with lower yields and 
fewer neutrons where Y(Af) is large to obtain the same V 
as the constant d^n where the neutrons from symmetric 
fission are effectively absent. 

Once the average total fragment kinetic energy has 
been determined, the average combined excitation en- 
ergy in the two fragments follows automatically by en- 
ergy conservation, 

Qlh - TKE = TXE = ~E* L + ~E* H ■ (15) 

The last relation indicates that the total excitation en- 
ergy is partitioned between the two fragments. The vari- 
ation of the total mean excitation energy with fragment 
mass is similar to that of v(A) in Fig. [9] 

FREYA assumes that the excitation energy is partitioned 
statistically, as it would be if the two fragments were 
in mutual thermal equilibrium. Consequently, TXE is 
divided in proportion to the heat capacities of the nascent 
fragments, 

E* = _ TXE, (16) 

ah + an 
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FIG. 8: (Color online) The average fragment kinetic energy 
as a function of fragment mass from Refs. [2?], H^] compared 
to FREYA calculations with a constant tip separation of do = 
4.05 fm and the average distance extracted from Fig. [7] 
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FIG. 9: (Color online) The average neutron multiplicity as a 
function of the fragment mass from Refs. [13, [II [If compared 
to FREYA calculations with a constant tip separation of do = 
4.05 fm and the average distance extracted from Fig. [7] 



where hi is the level-density parameter for fragment i. To 
take account of the microscopic structure of the individ- 
ual fragments as well as any possible energy dependence, 
FREYA uses the functional form due to Kawano et al. [3a| , 



~ ai (E*) = 

eo 



1 



Ui 1 



(17) 



where U% — E* — Aj and 7 = 0.05 |9|. The pairing energy 
of the fragment, A<, and its shell correction, SWi are 
tabulated in Ref. [Hj] based on the mass formula of Koura 
et al. [39J. Although FREYA uses the default value eo = 
7.25 MeV [2{|, we wish to make this value adjustable, 
taking 



e = (7.25 MeV) a 



(18) 



and treating the common factor a as a model parame- 
ter. We note that if the shell corrections are negligible, 



8W w 0, then this renormalization is immaterial and 
the excitation energy will be shared according to mass, 
E* oc Ai. 

The relationship between excitation energy E* and the 
temperature Tj is given by 



E* = cuTf 



(19) 



so that when the total excitation energy is shared ac- 
cording to the level-density parameters ai then the two 
fragment temperatures are equal, Tl = Th- 

While the equal temperature assumption is a reason- 
ably good first approximation, it may be inadequate for 
obtaining a detailed description of prompt neutron emis- 
sion. Therefore we redistribute the excitation energies of 
the fragments, 



E* L = xE* L 



E* H = TKE 



E 



L ■ 



(20) 



and treat x as an adjustable model parameter. The data 
indicate that the light fragments acquire more than their 
"fair share" of the energy, thus we expect that our sta- 
tistical analysis will favor x > 1. 

After the mean excitation energies have been assigned, 
FREYA considers the effect of thermal fluctuations in the 
partitioning of the excitation energy. For this task, FREYA 
assumes that the fluctuation in the excitation energy of 
a nucleus is a 2 E = 2E*T where T is its temperature and 
E* = aT 2 its mean excitation. Therefore, for each of the 
two fragments, we sample a thermal energy fluctuation 
SE* from a gaussian distribution of variance of = 2E*T{ 
and modify the fragment excitations accordingly, 



E* = E*+SE* , i = L,H 



(21) 



Due to energy conservation, there is a compensating op- 
posite fluctuation in the total kinetic energy, so 



TKE = TKE - SE* L - SE* H 



(22) 



With both the excitations and the kinetic energies of 
the two fragments fully determined, it is an elementary 
matter to calculate the magnitude of their momenta and 
thus sample the velocities with which they emerge after 
having been fully accelerated by their mutual Coulomb 
repulsion [29| . 



C. Neutron evaporation 

The primary fission fragments are typically sufficiently 
excited to permit the emission of one or more neutrons. 
For each of the two fragments, neutron emission is treated 
by iterating neutron evaporation from each fragment. 

At each step in the evaporation chain, the excited 
mother nucleus Ai Zi has a total mass equal to its ground- 
state mass plus its excitation energy, M* = Mf s + E* . 
The Q-value for neutron emission from the fragment is 
then Q n = M* — Mf — m n , where Mf is the ground-state 
mass of the daughter nucleus and m n is the mass of the 
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neutron (for neutron emission we have Af = Ai — 1 and 
Zf = Zi). The Q- value is equal to the maximum possible 
excitation energy of the daughter nucleus, which occurs 
if the final relative kinetic energy vanishes. The temper- 
ature in the daughter fragment is then maximal. Thus, 
once Q n is known, one may sample the kinetic energy of 
the evaporated neutron. FREYA assumes that the angular 
distribution is isotropic in the rest frame of the mother 
nucleus and uses a standard spectral shape [40| . 

f n (E) = =^ <* Ee-v/r- , (23) 
v dE 

which can be sampled very fast [29| . 

Although relativistic effects are very small, we take 
them into account in order to ensure exact conservation 
of energy and momentum, which is convenient for code 
verification purposes. We therefore take the sampled en- 
ergy E to represent the total kinetic energy in the rest 
frame of the mother nucleus, i.e. it is the kinetic energy of 
the emitted neutron plus the recoil energy of the residual 
daughter nucleus. The excitation energy in the daughter 
nucleus is then given by 

E* d = Q n - E . (24) 

The mass of the daughter nucleus is thus M d = Mf+E* v 
It is possible to calculate the magnitude of the momenta 
of the two final bodies: the excited daughter and the 
emitted neutron. Sampling the direction of their relative 
motion isotropically, we thus obtain the two final mo- 
menta which are subsequently boosted into the overall 
reference frame by the appropriate Lorentz transforma- 
tion. 

This procedure repeated until no further neutron emis- 
sion is energetically possible, when E* d < S a , where S n is 
the neutron separation energy for the daughter nucleus, 
S n = M( A *Z d ) - M{ A ^Z d ) - m n . 

V. RESULTS 

We now proceed to discuss our analysis. We first 
describe the computational approach and then explain 



how the model parameters are determined. The result- 
ing prompt neutron spectrum is then discussed in detail. 
Finally, we present some additional observables of par- 
ticular relevance. 



A. Computational approach 

FREYA is used to generate a large sample of fission 
events (typically one million events for each parameter 
set). For each set m of such randomly selected model 
parameter values, {s^ m \ a^ m \ x^}, the prompt fission 
neutron spectrum and v in each event m are then com- 
pared to the available experimental data at the given 
incident neutron energy, E n . This allows us to assign the 
likelihood for that particular set (see Sec. [TTJ) based on 
either the \m f° r comparison with 77 only, \% or on the 
total \m characterizing the comparison with both 77 and 
the spectral shape f n (E) = V^dV/dE, x| + Xspcctra: 



Wm = w{s {m \a {m \x^} = e- x ™/ 2 . (25) 



Since the weight w m depends exponentially on x^, the 
likelihood tends to be strongly peaked around the favored 
set. It is important that the parameter sample be suf- 
ficiently dense in the peak region to ensure that many 
sets have non- negligible weights. We typically sample 
2000 different parameter sets but have verified that the 
results remain unchanged when a five times larger sample 
is explored. 

Using this method, we can obtain those values of s, 
a and x that minimize either \% or Xv + Xspcctra- We 
denote the optimal set by {s°, a , x }. We also obtain the 
corresponding correlation matrix, as described in Sec. HU 



E n (MeV) 


S ° 


a 


x° 


V 


2 
Xu 




Xspcctra/-^ 


0.5 


1.05449 ± 0.00567 


1.10562 ±0.07987 


1.10264 ± 0.05909 


2.948 ±0.015 


4.26 x 10" 


-3 


28.99 


1.5 


1.05887 ± 0.00585 


1.10426 ±0.07854 


1.10178 ± 0.05736 


3.090 ±0.015 


8.46 x 10" 


-4 


9.81 


2.5 


1.06590 ± 0.00858 


1.10243 ± 0.07972 


1.09969 ±0.11359 


3.242 ±0.016 


1.88 x 10~ 


-2 


3.40 


3.5 


1.06886 ± 0.00902 


1.10440 ± 0.07903 


1.09987 ±0.11745 


3.373 ± 0.017 


3.78 x 10" 


-2 


5.90 


4.5 


1.07598 ± 0.00699 


1.10246 ± 0.07963 


1.09889 ± 0.05829 


3.527 ±0.017 


2.55 x 10" 


-2 




5.5 


1.08418 ± 0.00752 


1.10409 ± 0.08023 


1.09892 ± 0.05758 


3.681 ±0.019 


1.50 x 10~ 


-2 





TABLE III: The optimal values of the three model parameters s, a and x obtained in three-parameter fits to V alone, as well 
as the corresponding mean neutron multiplicities 77, together with the extracted uncertainties. The resulting values of Xtt an d 
Xspcctra per degree of freedom are also given. 



E n (MeV) 


s° 


a 


x° 


V 


xl 


Xspcctra 1 N 


0.5 


1.05705 ±0.00173 


0.96754 ± 0.02236 


1.00523 ± 0.00574 


2.961 ± 0.007 


0.76 


13.72 


1.5 


1.04573 ± 0.00742 


0.97291 ± 0.03424 


1.18356 ±0.05142 


3.078 ± 0.020 


0.43 


23.77 


2.5 


1.05485 ± 0.00602 


0.99909 ± 0.04221 


1.18587 ±0.06274 


3.239 ±0.016 


0.0066 


2.58 


3.5 


1.05309 ± 0.00657 


0.98038 ± 0.03839 


1.21052 ±0.05293 


3.364 ±0.013 


0.24 


4.61 



TABLE V: The optimal values s°, a 3 and x° obtained in three-parameter fits to the spectra and v. The corresponding values 
of V are also shown. The resulting \ 2 values for 77 and the spectra are given separately. 



B. Determination of the model parameters 

Tabic IIIII shows the optimal values and the associated 
uncertainties for the three model parameters used in our 
fission calculations. These values have been obtained by 
fitting only to the evaluated V while ignoring the spectral 
data. We have checked that fixing either x or a, or both, 
in these fits lead to equivalent results for all values of E n . 

The correlation coefficients between these model pa- 
rameters are shown in Table IIV1 If the parameters are 
uncorrelated, Ckk' = 0. Correlated parameters lead to 
nonzero correlation coefficients. If Ckk' > 0, a,k increases 
as cvfe' increases. On the other hand, if Ckk' < 0, a.k 
increases as ctk' decreases. The correlation coefficients 
between s and a, C sa , are relatively large and positive 
while those between s and x, C sx , are large and nega- 
tive, suggesting strong correlations between these pairs of 
parameters. The correlation coefficients between a and 
x, C ax , are close to zero and fluctuate in sign, signaling 
only a weak correlation between this pair of parameters. 
In contrast, when the spectra are also included in the 
fits for E n < 3.5 MeV, the correlation coefficients are all 
very close to ±1 in all cases, likely because the overlap in 
parameter space that simultaneously reproduces v and 
the spectra is small. 



E n (MeV) 


Cs a 


Cs X 


Ca x 


0.5 


0.608 


-0.569 


0.0156 


1.5 


0.611 


-0.561 


0.0042 


2.5 


0.465 


-0.776 


0.0212 


3.5 


0.464 


-0.766 


0.0441 


4.5 


0.757 


-0.569 


-0.0053 


5.5 


0.693 


-0.480 


-0.0130 



TABLE IV: The correlation coefficients (see Eq. (gj)) for the 
three parameters s, a and x fitted to V alone. 

The experimental values for the total average neutron 
multiplicity place remarkably stringent constraints on the 
value of the model parameter s while more room is left 
for variations of a and x. Specifically, changing the tip 



separation distance scale factor s by only 1% (keeping a 
and x fixed) changes v by 1.8%, far outside the experi- 
mental uncertainty. A change in s, see Eq. (|T4"|). results 
in a change in the average TKE, Eq. (fl"3")) . of less than 
0.5 MeV. Thus V is very sensitive to the balance between 
the kinetic and excitation energies. On the other hand, 77 
is less sensitive to the partition of the excitation energy 
between the light and heavy fragments since changing x 
by 5% (keeping a and s fixed) changes V by only 0.5%. 
Finally, 77 is least sensitive to changes in a which mod- 
ifies the fragment temperature, predominantly affecting 
the low energy part of the neutron spectrum. Changing 
a by 5% (keeping s and x fixed) changes 77 by only 0.3%. 

Table [V] shows results calculated by fitting to both 77 
and the prompt neutron spectra. (We do not show the 
4.5 and 5.5 MeV results again since there are no pub- 
lished spectra at these energies.) When the spectral data 
are included in the fit the agreement with these data 
and the evaluated v is poor. If we had confidence in the 
spectral data, this would be a formal indication that our 
model was incorrect or that uncertainties in 77 were un- 
derestimated. Inconsistencies in the spectral data (see 
Sect. IIIIBI) make either conclusion difficult. Some sets 
(particularly those of Alexandrova [20j, which make the 
largest contribution to the spectral \ 2 ) are inconsistent 
with other sets, and, in a number of cases, uncertainties 
conservatively estimated. In addition, the relative nor- 
malization, while determined from fitting to a Watt spec- 
trum and used only for scaling purposes, may increase 
the relative \ 2 for some data sets, possibly including the 
Aleksandrova sets which are only available for E > 1.5 
MeV. Indeed, since these sets give the largest contribu- 
tion to the total x 2 , eliminating them can change the 
optimal parameter values, while removing one or more 
of the other sets has little to no effect. For these rea- 
sons, we did not use the spectral data to obtain our fi- 
nal evaluation. In addition, as discussed in more detail 
later, there are indications from 235 TJ measurements that 
more neutrons are emitted from the light fragment than 
are from the heavy fragment (x > 1) [27j • The fit at 
E n —0.5 MeV shown in Table IVl is consistent with x = 1, 
giving V L KiV H - 



C. The prompt neutron spectrum in Fig. 1101 The top panel of this figure gives the 

A comparison between experimental data and our 
calculations of the prompt neutron spectrum is shown 



spectral shap e and shows all experimental data from 
Refs. [H OS OS M, HH H3, [H, M, [H . Since the shape 
varies slowly with incident neutron energy, the calcula- 
tions using parameters fit to V alone and to V and the 
spectral data are practically indistinguishable on a linear 
scale. The bottom panel of Fig. [TU] shows only the more 
recent Staples data from 0.5 to 3.5 MeV [2J|]. In this 
panel, the different spectra can be distinguished because 
they have been normalized to F, which varies modestly 
with incident neutron energy. 

Because FREYA cannot produce sufficient statistics at 
the fine energy scale needed by typical spectral evalua- 
tions, high statistics FREYA runs have been made to em- 
phasize the low and high energy tails of the spectra. To 
remove statistical noise, Watt distributions are fit to the 
low (E < 2 MeV) and high (E > 4 MeV) energy parts 
of the spectrum for each incident neutron energy. A fine 
grid is obtained in the intermediate part of the spectrum 
by interpolation. 

Figure [TT] gives the difference between the present cal- 
culations and the evaluations in ENDF/B-VII. Our spec- 
tra are systematically softer, giving lower mean neutron 
energies. This difference has important implications for 
criticality. 

In the previous section, we argued that currrently 
available spectral data should not be used in the fis- 
sion likelihood analysis. To illustrate the impact of these 
data on spectral calculations, we show the difference be- 
tween the fits without and with the spectral data at 
E n = 0.5 MeV, normalized to 77 on a log-log scale, in 
Fig. [T2J The difference is largest in the high-energy tail 
of the spectra where the fit to the spectra and v is softer. 
The ratio of the fits with and without the spectral data 
are shown in Fig. [131 Below 2 MeV, the fit with the spec- 
tral data is 1-2% higher but by E » 10 MeV, it is about 
60% lower than the spectral description with a fit to v 
alone. At higher energies the calculations grow further 
apart but the ratios are statistics limited since, even with 
1-2 million events FREYA does not fully populate the high 
energy tail of the emission spectrum. 

We can compute the uncertainty in the spectra as well 
as in the employed values of the model parameters. Each 
particular set of model parameter values, {a k }, yields a 
different neutron spectrum (dv/dE)^ so that the result- 
ing ensemble of spectra can be subjected to a statistical 
analysis at each value of the neutron energy E, yielding 
an average value of the neutron spectrum, dv/dE, and 
an associated dispersion, a(dvjdE). For E n — 0.5 MeV, 
Fig. HU shows the ratio between dvjdE + aidv jdE) and 
dv/dE. This spectral ratio provides an indication of 
the relative uncertainty on the spectrum at each energy. 
With the fits to 77 alone, the calculated uncertainty is 
less than 5% for E < 4 MeV and less than 2% for E < 2 
MeV, much smaller than the spread in the data depicted 
in Fig. 1101 The uncertainty increases approximately lin- 
early with E for E > 2.5 MeV, reaching « 40% at 15 
MeV. We have also shown the relative uncertainty with 
all spectra included in the fit as well as that obtained by 



11 




Neutron energy E (MeV) 




Neutron energy E (MeV) 



FIG. 10: (Color online) The measured prompt neutron spec- 
tra are compared to our fit results. The comparison to the low 
energy results from Refs. [13, El, E3, [13, [HI, [23, H3, [H, [ll] (tip- 
per panel) are of the normalized spectral shapes while the re- 
sults at higher incident neutron energies from Ref. }23j (lower 
panel) are compared to the spectral distributions themselves. 

leaving out the two spectra with the largest \ 2 . Both of 
these give small but noisy uncertainties, suggesting that 
result is not a true measure of the calculated uncertainty 
in this case and that the spectral uncertainty as shown is 
rather random. The noisiness of the combined fits is due 
to the difficulty of obtaining a combination of parameter 
values that simultaneously minimizes \% an d Xspectra- 

It is instructive to consider the correlations between 
the spectral strength at different energies. The eval- 
uation of the corresponding covariance (see Eq. ([6])) is 
complicated by the fact that the observables considered, 
specified energies of emitted neutrons, form a continuum. 
In practice, it is convenient to consider discrete energy 
bins (so the observable a k represents the mean number 
of neutrons emitted with a kinetic energy in the bin k 
centered at the energy value E k . Using Eq. (@|, we may 
then calculate the corresponding covariance matrix 

a(E k ,E k ,) = < (E k - E k ){E k , - E k ,) y . (26) 

However, it is important to recognize that for contin- 
uous observables, the above matrix function is singular 
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along the diagonal [42| . 

a(E kl E k ,) = 5| 4 <5(£ fe - E k ,) + a Ek E k , , (27) 

where a\ is the variance in the differential yield at the 
specified energy E k , while &E k E k , expresses the correla- 
tion between the differential yields at two different en- 
ergies E k and E k >. To obtain this latter quantity, we 
must first remove the singular part. This can be readily 
accomplished when the observable has been discretized 
by simply replacing the diagonal elements in a(E k} E k r) 
by values obtained by interpolating between the near- 
diagonal elements. The resulting correlation coefficient, 

C(E k ,E k ,) = o Ek Ej[°E k °E k \ , (28) 

is then regular. It is displayed in Fig. [T5] for the ensem- 
ble obtained for E n — 0.5 MeV by fitting to V alone. 
Figure [16] shows cuts at constant total neutron energy, 
E k + E k i . Similar results are found for all other incident 
energies considered. 

When the model parameters are varied, the spectral 
shapes tend to pivot around £si2 MeV. Consequently, 
when both neutron energies lie on the same side of this 
value, the differential changes are in phase and the cor- 
relation coefficient is close to one.O The changes are in 
opposite directions when the two energy values are on 
opposite sides of the pivot energy. By contrast, when the 
spectral data are included in the fits, the correlation co- 
efficients vary widely between +1 and —1 in no apparent 
pattern. 

When the number of FREYA events included in the 77- 
only fits at E n = 0.5 MeV is increased by a factor of 
five, the fitted model parameter values change by less 
than one standard deviation. When the spectra are also 
included in the fits, the resulting change in the fitted pa- 
rameter values increases Xv from 0.75 to ~15 without 
significantly improving the spectral fits. Moreover, while 
the fluctuations in the energy correlation coefficients de- 
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FIG. 12: (Color online) Prompt neutron spectra calculated in 
the laboratory frame as a function of outgoing neutron energy 
for 0.5 MeV incident neutron energies. The solid curve is 
obtained by fitting V alone while the dashed curve is fit to 
both the spectra and V. 
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FIG. 13: (Color online) The ratio of the spectra obtained by 
fitting to V and the spectral data relative to a fit based on V 
alone at E n = 0.5 MeV. 
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FIG. 11: (Color online) The percent difference between our 
evaluated spectra and that of ENDF-B/VII for all six incident 
neutron energies considered. 



crease somewhat when the larger event samples are used, 
they do not disappear. 

As is the case for the model parameters, there are un- 
certainties in the spectral calculations. If the model is 
qualitatively wrong, and the right spectral form cannot 
be obtained by simply changing the parameter values, 
then the spectral uncertainties are not correct. To ex- 
plore this we performed several additional variations on 
the model. In Sec. HVI we saw that a model that employs 
a constant tip separation d, independent of the specific 
binary partition, reproduces neither the total kinetic en- 
ergy data nor the neutron yield as a function of fragment 
mass. Nevertheless, a constant d^n yields better agree- 
ment with the 77-only fit in Fig. [T^] than with the fit that 
also includes the spectra. Similarly, making the level den- 
sity parameter independent of energy, a — A/e , changes 
the spectrum by less than 5% at lower energies (< 5 MeV) 
and by less than 20% at higher energies. Fundamental 
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FIG. 14: (Color online) The spectral ratios (see text) for the 
three different analyses indicated. 
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Neutron energy E1 

FIG. 15: (Color online) The correlation coefficients, 
C(Ei, E2), for the spectral strength of the evaporated neu- 
trons (see Eq. (|28[0 . as obtained from the statistical analysis 
at E„ — 0.5 MeV when only V is considered in the fits. Fig- 
ure [16] shows cuts along the three indicated lines of constant 
total energy. 

microscopic calculations of fission could provide insight 
into the sensitivity of the spectrum to changes in the 
parameters, leading to better estimates of the spectral 
uncertainties. 

Critical assemblies, which are designed to determine 
the conditions under which a fission chain reaction is sta- 
tionary, provide an important quality check on the spec- 
tral evaluations. The key measure of a critical assembly 
is the neutron multiplication factor fc e gf (often denoted 
as the k eigenvalue). When this quantity is unity, the 
assembly is exactly critical, i.e. the net number of neu- 
trons resulting from each neutron-induced fission event 
is one on the average. (This number is the difference be- 



tween the number of neutrons emitted during the fission 
process and those lost to absorption and escape.) The 
degree of criticality of a particular assembly depends on 
the multiplicity of prompt neutrons, their spectral shape, 
and the (n, f) induced-fission cross section. 

Plutonium criticality is especially sensitive to the 
prompt neutron spectrum because the 239 Pu(n, /) cross 
section rises sharply between E n — 1.5 and 2 MeV. As a 
result, increasing the relative number of low-energy neu- 
trons tends to decrease criticality, lowering fc e ff, while 
increasing the number of neutrons having higher energy 
increases criticality. 

Figure [T7] shows calculations of fc e ff for different pluto- 
nium assemblies. Apart from the spectra, all data used in 
these calculations were taken from ENDF/B-VII. Overall 
there is good agreement with the measured values of fc e ff, 
though this new softer spectrum decreases the calculated 
values by about 0.003. Since this is approximately 1.5 
standard deviations away from the measurement, there 
may be an indication that the Pu fission cross section 
or neutron multiplicity is low by about a tenth of a per- 
cent. There appears to be room for some adjustment of 
the experimental data since the uncertainties in the cross 
sections are about 1%, while those in V are about 0.5%. 



D. More exclusive observables 

Though less important for understanding energy pro- 
duction, more exclusive observables play a central role 
in the development of a comprehensive description of fis- 
sion. Figure [18] shows calculations of fragment kinetic 
and excitation energies. Note that the fragment kinetic 
energies are almost independent of the incident neutron 
energy. Indeed, the kinetic energy appears to decrease 
slightly with energy, as may be expected since s increases. 
This may at first appear surprising but the Coulomb ap- 
proximation to the total kinetic energy in Eq. (fT3| is 
independent of the incident neutron energy. These re- 
sults are consistent with measurements made with 235 U 
and 238 U targets over a similar incident neutron energy 
range, 0.5 < E n < 6 MeV O and 1.2 < E n < 5.8 MeV 
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FIG. 16: (Color online) The spectral correlation coefficients, 
C(Ei, E2), along the three lines of constant combined energy 
indicated in Fig. \T5\ 
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FIG. 17: (Color online) Calculated k e g for several 239 Pu crit- 
ical assemblies obtained using our fits for 0.5 < E„ < 5.5 
MeV in the Mercury Monte Carlo. The results are compared 
to those with the standard ENDL2008.2 and ENDF-B/VII 
databases. 



[44[ respectively. In both cases, the average TKE, TKE, 
changes less than 1 MeV over the entire energy range. 
Ref. [4J| also shows that, while the mass-averaged TKE 
is consistent with near energy independence, higher en- 
ergy incident neutrons typically give more TKE to masses 
close to symmetric fission and somewhat less TKE for 
Ah > 140. The slight increase in TKE close to symmet- 
ric fission of 238 TJ is not unexpected since the symmetric 
contribution to Y(A) increases with incident neutron en- 
ergy. Since such detailed TKE information is not avail- 
able for neutrons on 239 Pu, we have therefore chosen to 
use a constant scale factor at each energy. 

The constancy of the fragment kinetic energy as a func- 
tion of E n allows the energy of the incident neutron to be 
converted into excitation energy. The increase of E* with 
E n is fairly monotonic over all Af , see the right-hand side 
of Fig. [TH It appears, however, that the slope of v(Af) 
for Af > 120 increases somewhat faster with E n than 
for Af < 120, as shown on the left-hand side of Fig. \T§\ 
This scenario is consistent with washing out the saw- 
tooth pattern of v(Af) with increasing neutron energy 
|45| . See Table I VII for the average neutron multiplicity 
for the light and heavy fragments as well as the sum. The 
associated multiplicity dispersions, a v = [(v) 2 — V 2 ] 1 / 2 , 
are also given. Since V is used to determine the values 
model parameters, it may be preferable to use a different 
(and more exclusive) observable to check whether a given 
model parameter set is preferred over another. A better 
choice is the average neutron multiplicity and average 
neutron energies from the individual fragments. There 
are some limited data on thermal neutron-induced fis- 
sion of 235 U and spontaneous fission of 252 Cf [47| 
which suggest that the light fragment emits more neu- 
trons than the heavy fra gme nt, 40% more for 235 U [46| 
and 20% more for 252 Cf 0. Our results for 0.5 MeV, 
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FIG. 18: (Color online) The average kinetic energy of the 
fission fragments {upper panel) and their average excitation 
(lower panel) as a function of fragment mass number Af for 
0.5 < E n < 5.5 MeV. 



shown in Table fVTl give a relative difference in 17 between 
the light and heavy fragments of about 20% for x ~ 1.1. 
Fits to v and the spectral data rather than v alone give 
vl ~ vh for E n < 0.5 MeV, seemingly excluded by these 
measurements, if the same is true for Pu. 

A more sensitive neutron observable is the kinetic en- 
ergy of an evaporated neutron. The lower panel in Fig. 1191 
shows the average kinetic energies of the emitted neu- 
trons as a function of fragment mass for the lowest and 
highest incident energies studied (0.5 and 5.5 MeV). The 
average kinetic energy of the emitted neutrons is almost 
constant with A except in the region 110 < A < 140 
where it increases. The dip in TKE occurs in the sym- 
metric region, making more energy available for neutron 
emission, resulting in more and faster prompt neutrons. 

In Fig. [20] we show the probability for a given neutron 
multiplicity, P(v), as a function of neutron number for all 
E n . Along with the probability distribution for emission 
from both fragments, we also show the distributions for 
the light and heavy fragments separately. 

Table IVIII gives the average energies of the neutrons 
emitted from the light fragment, the heavy fragment, or 
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FIG. 20: (Color online) The normalized neutron multiplicity 
distribution obtained with FREYA for both fragments (top), 
the light fragment (middle) and the heavy fragment (bottom) . 
Results are shown for 0.5 < E n < 5.5 MeV. 
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FIG. 19: (Color online) The average neutron multiplicity for 
0.5 < E n < 5.5 MeV (upper panel) and the mean kinetic 
energy of the evaporated neutrons (lower panel) for E n = 0.5 
and 5.5 MeV, as functions of the fission fragment mass number 



from either, together with the associated variances, for 
the incident neutron energies E n . The average energies 
increase with E n in all cases and those coming from the 
light fragment tend to be more energetic than those com- 
ing from the heavy one, so we have (E L ) > (E L+H ) > 
(E H ). The variances exhibit the same hierarchy as the 
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TABLE VI: The mean combined neutron multiplicities V as 
well as the mean multiplicities of neutrons emitted from either 
the light or the heavy fragment, Vl and Vh, together with the 
associated dispersions. 



TABLE VII: The average energy of neutrons emitted by either 
fragment and by the light and heavy fragments separately, 
along with the associated dispersions, for various incoming 
neutron energies (all in MeV). 



average energies but increase more slowly with incident 
energy. The overall average energy (E) is similar to that 
obtained for thermal neutron-induced fission of 235 U and 
252 Cf(sf) found in Ref. 0]. 



VI. CONCLUSION 

Our studies employ the recently developed a Monte- 
Carlo model, FREYA, that simulates fission and the sub- 
sequent neutron and photon emission from the fragments 
on an event-by-event basis, maintaining energy and mo- 
mentum conservation at each step in the production and 
de-excitation of the fragments. We have introduced three 
adjustable parameters, s, a, and x, which modulate the 
separation between the tips of the fragments, scale the 
level-density parameter for the fragments, and modify 
the partition of energy between them, respectively. These 
three model parameters were varied over an appropriate 
range and, for each particular set of values, FREYA was 
used to generate a large sample of fission events from 
which the resulting properties of the neutron spectra were 
extracted. Each set of parameter values was assigned a 
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likelihood weight based on the x 2 obtained from compar- 
ison with the measured mean multiplicity V and/or the 
measured differential neutron spectrum dv/dE. Mean 
values and covariances for both input parameters and 
quantities predicted by the model were obtained through 
standard statistical techniques. This combination of the 
Monte-Carlo fission model with the likelihood weight- 
ing presents a powerful tool for the evaluation of fission- 
neutron data. 

This procedure was applied to the analysis of neutron- 
emission data for neutron-induced fission on 239 Pu, from 
thermal to 5.5 MeV incident energies. Although the ap- 
proach taken and the nucleus studied in this work are 
different, the results largely corroborate the findings of 
Lemaire et al. Q in emphasizing the importance of the 
initial conditions {e.g. the kinetic and excitation energies 
of the fragments). Furthermore, our work underscores 
the effectiveness of the measured V in constraining the 
model parameters, more strongly even than the differen- 
tial neutron-spectrum data. In particular, it was found 
the the parameter controlling the tip separation between 
fragments was by far the most important in reproducing 
the experimental V values. In the end, fits of our model to 
the v data alone (i.e., excluding the differential-spectral 
data) were found to be more robust and were used to 



obtain the best model parameters. 

We plan to apply this method to the prediction of 
neutron emission properties in other actinides. However, 
in those cases where critical experimental data (such as 
kinetic energies and neutron multiplicities and spectra) 
are not available to constrain the FREYA calculations, it 
may be necessary to invoke supplementary information 
from various theoretical models, such as Hartree-Fock- 
Bogoliubov or macroscopic-microscopic treatments. 
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